跟Nature一起学习TCGA,GTEx和CCLE数据库的使用
这是果子学生信文献阅读群的第一个产品，起源可以看一下这篇帖子，
科研工作者的力量来自持续不断的读文献！
这一篇文献是陈建军发表在Nature上的letter
Histone H3 trimethylation at lysine 36 guides m6A RNA modification co-transcriptionally
这篇Nature首次揭示了组蛋白修饰对于m6A甲基化的影响，展示了基因表达调控的新模式。但是这篇文章让我念念不忘的是文中的一张附图 。
因为SETD2是H3K36me2到H3K36me3特异的甲基转移酶，而MTC复合物（包含METTL3, METTL14和WTAP）又是m6A甲基化的转移酶。如果H3K36me3和m6A有关系，那么SETD2和MTC成员也有可能表达上相关。在没有做实验之前，作者是通过使用3个大型公共数据库的数据来说明的。
这三个数据库分别是GTEx，TCGA，CCLE。
其中GTEx数据库中是正常人体各个组织的转录组数据
TCGA是美国的癌症计划，里面包含了33种癌症组织的多组学数据
CCLE中有1000多个癌症细胞系，都做了转录组测序。
作者使用这三个数据库来说明SETD2和MTC的三个成员表达正相关，其中CCLE中的图，就是图e，是我们能理解的相关性分析，横坐标一个基因，纵坐标一个基因，打点然后画出拟合线。但是GTEx和TCGA中的相关性图我没见过。
以图d蓝色圈线中的图为例，横坐标是p值取log后的负值，这样越往右边越大，竖线应该是0.05的位置，横坐标是相关性系数，横线是0。我立马就知道，这是多个癌症的相关性结果的汇总图。因为TCGA有33个肿瘤，每个肿瘤都很很多样本，两个基因在单个肿瘤中计算相关性会得到P值和相关性系数，把33个癌症中的结果汇总就得到了图d，这个本质上就是个火山图。
这个图可以帮我们很好地展示普适性结果，比如，图d可以看出SETD2和METTL3在大部分癌症中都是正相关的，在个别肿瘤中不相关。同理GTEx因为记录了人体多个组织的信息，也可以画出类似的图。
当我看到这个图的时候，很高兴，因为这是张科研人员人人都用得到的图。医学科学研究中，大部分情况都是关注两个分子的相关作用，尤其是转录调控，如果转录因子A调控了靶基因B, 那么A和B就应该是个正相关的关系，当我们有多个转录因子可以选择的时候，用这里的方法可以缩小范围。
因为这个生信技能太实用，我花了2天时间（其中一天是下载数据）像素级别还原了这里的9张图。而且，方法具有普适性，给出任意两个基因，都可以很方便的出图，人人可用！
以下是教程的思路，也是我整理数据的常用流程，熟悉R语言后很快就会完成，R语言不熟悉的文末有视频教程。这里主要是三个步骤
第一，下载数据
第二，转换数据: 清洁数据三部曲，基因注释，行列转置，添加分组
第三，计算作图
我们以TCGA的数据为例。
首先是下载数据
https://xenabrowser.net/datapages/?cohort=TCGA%20Pan-Cancer%20(PANCAN)&removeHub=https%3A%2F%2Fxena.treehouse.gi.ucsc.edu%3A443
打开链接找到这个部分
点进去下载表达量数据
再下载临床数据作为分组
下载gtf注释文件，该数据使用的是gencode23版本
https://www.gencodegenes.org/human/release_23.html
然后是清理数据
第一，读入数据
第二，用gtf文件转换基因名称
第三，用
t()行列转置第四，从临床数据里面找到分组信息并添加
这里面的subtype是TCGA id的第13，14位，我在这里去除了癌旁组织
这里type就是癌种信息，这里有33个，比如乳腺癌就用“BRCA”来表示
现在数据处理完了，
最后计算作图
这里我写了一个函数，批量计算各个癌症中的相关性系数和p值
pancor <- function(gene1,gene2,data){
data1 <- split(data,data$type)
do.call(rbind,lapply(data1, function(x){
dd <- cor.test(as.numeric(x[,gene1]),as.numeric(x[,gene2]),type="pearson")
data.frame(type=x$type[1],cor=dd$estimate,p.value=dd$p.value )
}))
}
以METTL3,SETD2这两个基因为例调用函数
plotdf <- pancor("METTL3","SETD2",data)
最终得到一个数据框，有三列，第一列是癌种简称，第二列是相关性系数，第三列是p值
现在开始作图
library(ggplot2)
options(scipen = 2)
ggplot(plotdf,aes(-log10(p.value),cor))+
## 画点
geom_point(size=6,fill="purple",shape = 21,colour="black",stroke = 2)+
## 改变y轴
scale_y_continuous(expand = c(0,0),limits = c(-1,1),breaks = seq(-1,1,0.2))+
## 改变x轴
scale_x_log10(expand = c(0,0),limits = c(0.1, 1000),breaks = c(0.1,10,1000))+
## 加横线
geom_hline(yintercept = 0,size=1.5)+
## 加竖线
geom_vline(xintercept = -log10(0.5),size=1.5)+
## 坐标轴
labs(x=bquote(-log[10]~italic("P")),y="Pearson correlation (r)")+
theme(axis.title=element_text(size=20),
axis.text = element_text(face = "bold",size = 16),
axis.ticks.length=unit(.4, "cm"),
axis.ticks = element_line(colour = "black", size = 1),
panel.background = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.line = element_line(colour = "black"),
panel.border = element_rect(colour = "black", fill=NA, size=1.5),
plot.margin = margin(1, 1, 1, 1, "cm"))
出的图是这样的
对比一下原图，有一点点的差别，这应该跟样本的选取数量还有数据的预处理有关，可以不用管。我们这里下载的是TPM的数据，而且是已经log2转换后的。
此时，我们可有实现任意两个基因的相关性作图,再换一个试试。
plotdf <- pancor("METTL14","SETD2",data)
因为这里的图是用ggplot2画的，可以根据自己的喜好任意改颜色和风格。
如果想要关注某个特定的癌症种两个基因的相关性，那就更简单了，把对应癌症的行选取出来，直接计算
这是乳腺癌中METTL3,SETD2这两个基因的相关性图，也是高度可定制的。
好了，GTEx和CCLE的数据也是同样的处理思路。R语言的老手搞定这些肯定都不在话下，但是对于无法独自搞定的朋友，我们给出了以下方案：
1.如果你有这样的需求，给我两个基因，我帮你画三张图，收费200RMB/次。
2.如果你想自己学习(推荐)，实现任意两个基因画图，我也录制了一个详细的视频教程，包含所有的输入文件和代码(R 语言project格式交付)，大概1个多小时，定价129元。
因为这是第一个产品，我设置了一个三天拼团的活动，活动价格为99元。（在微信的微店里面）
购买后请加我微信guotosky，我来拉个群答疑。
3.如果你没有接触过R语言，也想实现这个人人可用的操作，需要先完成一件事情。在果子学生信公众号回复“果子学生信”，有一个R语言安装和环境配置的教程，被大概1000个人检验过，跟着视频无脑操作即可。
4.如果你在安装的过程中出现无法自行解决的问题，我们还准备了一个线上远程服务，由我们的专业人员连线帮你搞定。收费200 RMB。
科研工作者的力量来自持续不断的读文献，没毛病。
阅读原文也有链接。